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We present details of numerical simulations of the gravitational radiation produced by a first 
order thermal phase transition in the early universe. We confirm that the dominant source of 
gravitational waves is sound waves generated by the expanding bubbles of the low-temperature 
phase. We demonstrate that the sound waves have a power spectrum with a power-law form between 
the scales set by the average bubble separation (which sets the length scale of the fluid flow Lf) and 
the bubble wall width. The sound waves generate gravitational waves whose power spectrum also has 
a power-law form, at a rate proportional to Lf and the square of the fluid kinetic energy density. We 
identify a dimensionless parameter Slew characterising the efficiency of this “acoustic” gravitational 
wave production whose value is SyrDow — 0.8 ±0.1 across all our simulations. We compare the 
acoustic gravitational waves with the standard prediction from the envelope approximation. Not 
only is the power spectrum steeper (apart from an initial transient) but the gravitational wave 
energy density is generically larger by the ratio of the Hubble time to the phase transition duration, 
which can be 2 orders of magnitude or more in a typical first order electroweak phase transition. 

PACS numbers: 64.60.Q-, 47.75.±f, 95.30.Lz 


I. INTRODUCTION 

Gravitational waves (GWs) promise a new and exciting 
window to the cosmos. Two ground-based interferometer 
experiments, LIGO and VIRGO, are about to restart op¬ 
erations with greatly increased sensitivity 0,0 , and will 
be joined in a few years by KAGRA 0. Once working at 
their design sensitivity, they are expected to quickly de¬ 
tect gravitational wave signals from binary neutron stars 
0 . Gravitational waves also offer a unique way to learn 
about the early universe. A range of phenomena, such as 
inflation, topological defects and phase transitions may 
lead to observable gravitational wave signals across a 
wide range of frequencies (for a review see 0). There 
are a number of proposals to realise a gravitational wave 
detector in space, in the first place eLISA 0, which is 
scheduled for launch in 2034. Space-based detectors have 
much longer arm lengths than ground based ones, and 
have maximum sensitivity in a frequency range which 
is relevant for a first order phase transition at the elec¬ 
troweak scale. 

Given these exciting observational prospects, we re¬ 
visit the generation of gravitational waves in first or¬ 
der thermal phase transitions in the early universe. We 
have in mind an electroweak-scale phase transition, but 
nothing in our formalism is specific to electroweak scale 
physics. In the Standard Model the electroweak transi¬ 
tion is known to be a cross-over 0E1 , which does not 
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lead to a gravitational wave signal. However, a strong 
first order phase transition is possible in various exten¬ 
sions of the Standard Model [l2l4l^ . 

We reduce the original physical system to a model con¬ 
sisting of a scalar order parameter field coupled to an 
ideal fluid. The parameters of the model can in principle 
be fixed by matching to the thermodynamical quanti¬ 
ties of the original theory. We perform very large scale 
numerical simulations to determine the fluid and gravita¬ 
tional wave power spectra. The ultimate goal is to under¬ 
stand what information on the phase transition can be 
extracted from the future observation of a gravitational 
wave signal. 

Since the early nineties there have been a number of 
studies of gravitational waves from phase transitions. In 
Refs. |19l - l^ . the case of a scalar field only, i.e. a vacuum 
transition without fluid, was considered, motivated by 
models of inflation terminated by a first order transition. 
Vacuum transitions during inflation and with a fluid were 
considered in Ref. [2^. 

In a vacuum transition, all the energy released goes 
into the bubble wall, which as a result is accelerated to 
the speed of light. After solving numerically the field 
equations for the collision of two scalar bubbles it 

was realised that the energy-momentum tensor sourcing 
the gravitational wave production can be app roximated 
by the envelope of the colliding bubbles [^, [^. This 
“envelope approximation” models a configuration of ex¬ 
panding bubbles by the overlap of a corresponding set 
of infinitely thin shells. The envelope disappears once 
the transition is completed and gravitational wave pro¬ 
duction stops. It is found that the gravitational wave 
spectrum peaks at a frequency determined by the aver¬ 
age bubble size at collision. In the UV, the spectrum falls 
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as a power law, subsequently shown to be k~^ 1^. w here 
k is the wave number. Numerical studies in Ref. [26l| did 
not have the dynamic range to clearly confirm this be¬ 
haviour, but the larger simulations done for the present 
work show some supporting evidence. 

The case of a thermal phase transition, where the 
scalar field is coupled to a fluid, is more complicated. The 
nucleated bubbles will show accelerated expansion until 
the pressure inside is balanced by friction caused by the 
plasma. The bubbles then expand with a constant veloc¬ 
ity. This is because the energy released by the transition 
grows with the volume of the bubble, i.e. ^ R^, while the 
energy transferred to the scalar bubble wall only grows 
with the bubble surface, i.e. ~ R^, where R denotes the 
bubble radius. Hence only a tiny fraction of the released 
energy, on the order of the ratio of the initial to final 
bubble radius, stays in the scalar field. In the case of 
a first order electroweak scale thermal phase transition 
this ratio is about ^ 10“^^. Therefore 

gravitational wave production in thermal phase transi¬ 
tions is completely dominated by the fluid HI The energy 
which is released into the fluid mostly goes into reheating 
the plasma. A small and calculable fraction [2^, [2^ goes 
into bulk motion of the fluid and can source gravitational 
waves. 

Having established the fluid as the main source of grav¬ 
itational waves, the question of the production mecha¬ 
nism arises. Several mechanisms have been suggested 
and studied in the literature. In the simplest approach, 
one assumes that the fluid put into motion by the scalar 
wall can still be treated as a thin shell and the energy mo¬ 
mentum tensor sourcing gravitational wave production 
can again be approximated by the shell overlaps [2^ . In 
this case gravitational wave production finishes with the 
completion of the phase transition, and a characteristic 
prediction is the k~^ UV power law of the spectrum. An¬ 
other possibility is that the collision of bubbles induces 
turbulent motion of the fluid [2^. The resulting eddies 
generate gravitational waves even after the transition is 
completed [H, . Various UV power laws of the the 

gravitational wave spectrum have been suggested in this 
context, such as k~^-^ [s^ and [s^. 

To shed light onto these competing scenarios, we re¬ 
cently performed large scale numerical simulations of a 
thermal phase transition of a scalar field plus fluid sys¬ 
tem [s^. We found no indications that fluid turbulence 
was an important source of gravitational radiation. In¬ 
stead sound waves are generated by the explosive bub¬ 
ble growth, which propagate through the plasma until 
long after the transition is completed. In our simulations 
these sound waves are the dominant source of gravita- 


^ An exception may be the case where the bubble wall “runs away”, 
i.e. friction is not sufficient to prevent the wall from approaching 
the speed of light [^ . similar to a vacuum transition. Then 
both the scalar and the fluid could contribute sizeably to the 
generated gravitational wave signal. 


tional waves. After the phase transition, the fluid energy- 
momentum tensor clearly does not show the form as¬ 
sumed in the envelope approximation. The nearly linear 
behaviour of sound waves is very different to the highly 
nonlinear behaviour of the scalar field. 

Other numerical simulations of the generation of grav¬ 
itational waves by the coupled field-fluid system, using 
an explicit update algorithm for the fluid, have been de¬ 
scribed recently in Ref. [s^. The generation of gravi¬ 
tational waves through sound in QCD and electroweak 
phase transitions was also recently studied in Ref. [s^ . 
with special focus on the effect of possible non-linear 
sound dispersion relations, which were argued to lead to 
an inverse acoustic cascade. In Ref. [13, generation of 
gravitational waves in the Standard Model in the absence 
of such a cascade was discussed. 

In the present work we simulate at larger volumes and 
larger average bubble separations than in [13 , for the 
same range of bubble wall speeds and phase transition 
strengths. We widen the dynamic range even more by 
nucleating all bubbles at the same time. We confirm that 
the gravitational wave density parameter is proportional 
to the fourth power of the mean square fluid velocity, the 
ratio of lifetime of the source to the Hubble time, and the 
ratio of length scale of the source to the Hubble length. 
We measure the length scale of the source, approximately 
the average bubble separation in [13, directly from the 
fluid flow. With this improvement, the proportionality 
constant for the gravitational wave density parameter 
varies much less between phase transitions with differ¬ 
ent strengths and bubble wall speeds. Our measurements 
show that it is 0.8±0.1, where the uncertainly is the root 
mean square fluctuation between simulations. 

We show that the resulting gravitational wave spec¬ 
trum exhibits UV power laws which are clearly steeper 
than the k~^ predicted by the envelope approximation. 
In the case of deflagrations (where the bubble walls are 
subsonic), we are reasonably confident that the power 
law is k~^. For detonations we do not have sufficient 
dynamic range to be certain of the power law index. 

We compare the acoustic gravitational waves with the 
standard prediction from the envelope approximation. 
We argue that the envelope approximation is based on an 
incorrect picture of the dynamics of the fluid, in which 
the fluid perturbations are destroyed by bubble collisions 
in the same way as the bubble walls. Instead, they pass 
through one another, and keep oscillating, resulting in a 
gravitational wave source whose effective lifetime is the 
Hubble time. The true gravitational wave energy density 
is therefore a factor /3/i?* higher, where if* is the Hub¬ 
ble rate at the phase transition, and is the duration 
of the phase transition. For a thermal electroweak-scale 
phase transition, the gravitational wave signal is larger 
than hitherto believed by at least two orders of magni¬ 
tude. 



3 


II. FIRST ORDER PHASE TRANSITIONS IN 
COSMOLOGY 


A. Hydrodynamics 


We describe the phase transition using the cosmic fluid 
- order parameter field model [s^, which we sum¬ 
marise here. The model contains a classical scalar field (j) 
(effective order parameter), which is coupled to ideal fluid 
hydrodynamics. The variables describing the local state 
of the matter are local temperature T, fluid 4-velocity 
and the scalar order parameter field (j). The first or¬ 
der dynamics are obtained by introducing a temperature 
dependent effective potential V{(j),T). Following 
we use a simple cj)^ form for the potential: 

T) = - ^-AT4,^ + (1) 

The detailed form of the potential is not important, 
as long as it allows for a first order phase transition 
of sufficient strength. A first order transition occurs if 
2A2 < 9Ay. 

The equation of state of the coupled scalar field and 
fluid system is 

r)V 

e(r,0) = 3ar4 + F(0,T)-T—, (2) 

p{T,cj,)=aT^-V{cj,,T) (3) 

where a = (7r^/90)g*, and g* is the effective number of 
degrees of freedom. The latent heat density (usually just 
called the latent heat) is 

C{T) = wiT,0)-w{T,M (4) 


where w = e+p is the enthalpy density, and (j)ti is the equi¬ 
librium value of the field in the symmetry-broken phase 
at temperature T. The strength of the transition can be 
characterised by the ratio of the latent heat to the to¬ 
tal radiation density in the high temperature symmetric 
phase. 


^ 3ar4 ■ 


( 5 ) 


The total energy-momentum tensor of the system can be 
written as 


- \g^''{d4>f + [e+p] + g^^’^p (6) 

where the metric convention is (-H—h-b). The energy- 
momentum tensor is conserved, = 0. The interac¬ 

tion between field gradients and the fluid is introduced 
by splitting the conserved current nonuniquely into field 
and fluid parts, which are then coupled together through 
a dissipative term proportional to field gradient: 

= S'' (7) 

= 5^[(e + pW^^U''] - d-'p + 

( 8 ) 


where the coupling term is 

S'' = (9) 

with rj an adjustable friction parameter [s^. Equations 
analogous to Eqs. Q-® can, at least in principle, be de¬ 
rived from field theory (see e.g. [dll,[d^), but the simpli¬ 
fied model here is adequate for parametrising the entropy 
production (d^ . 

From Eqs. o and ® we can derive the equations of 
motion in a form suitable for numerical simulation. For 
the field we obtain 

dV 

+ —= 7?W(</> + W5,<(.), (10) 

where W is the relativistic y-factor and V is the fluid 
3-velocity, W = WV. For the fluid energy density E = 
We, contracting with C/^ gives 


E + d,{EV^)+p[W + d,{WV^)] + 

= gW^{^ + V^d,cj^f. ( 11 ) 


Finally, the equations of motion for the fluid momentum 
density Zi =W{e +p)Ui are 


Z, + dj{ZW^) + d,p + —5,0 = -gW{<p + 0)5,0. 

( 12 ) 

The implementation of Eqs. (Iini)-(II2|) on a discrete lattice 
is described in Section IlYl 

The parameters of the potential in Eq. o are related 
to thermodynamic quantities at the phase transition: the 
critical temperature Tc, latent heat C{Tc), surface tension 
a and the broken phase correlation length (which is also 
of order the bubble wall thickness) i [4^ 


C 

a 

e 


-^0 

1 - 2AV(9Ay) 

(13) 

rj,2rr2 

A2 

(14) 

2 V 2 A3 

81 A5/2 ^ 

(15) 

9A 1 

(16) 

2A2 ■ 


Due to supercooling, the phase transition (bubble nucle- 
ation) starts at temperature T^r, where Tq < < Tc. 

We are mostly interested in the large supercooling (LSC) 
case, where T/v is typically somewhere in the middle be¬ 
tween Tq and Tc- However, we emphasise that our fo¬ 
cus in this work is not the nucleation of critical bubbles, 
which in a given microscopic theory is a thermal field the¬ 
ory problem and can be studied in perturbation theory 
or with numerical simulations 0. In our simulations 
both the density of the initial bubbles and the nucleation 
temperature Tn are set by hand. 
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B. Bubble nucleation 


III. THEORY OF GW GENERATION 


The phase transition proceeds by the nucleation and 
growth of bubbles of the broken phase [1^, . Bubble 

nucleation occurs at an exponentially grow ing rate per 
unit volume below the critical temperature [ifl , 


p{t) ~ 


(17) 


A. GW power spectrum definition 

A gravitational wave is a propagating mode of the 
transverse and traceless part of the metric perturbation, 
hij. We are interested in calculating the gravitational 
wave energy density power spectrum, where the gravita¬ 
tional wave energy-momentum tensor is 


where —/? is the time derivative of the action of the criti¬ 
cal bubble S{t), and Tq is a dimensional prefactor of order 
c^w'^c Hli where aw ~ 1/30. The nucleation time In 
can be defined to be the time at which the nucleation 
rate reaches one bubble per Hubble volume per Hubble 
time, or p{t^) = 

The tunnelling rate parameter /3 not only sets the 
timescale of the transition, but also the average separa¬ 
tion between bubbles once the transition has completed, 
i?*. Having defined i?* to be the inverse cube root of the 
number density of bubbles, it can be shown that [1^ 

i?. = (87r)^^. (18) 


^ {d,h,Ah.o) ■ ( 21 ) 

To this end, we define the spectral density of the time 
derivative of the metric perturbation P^(k,f) by 

(/iy(k,<)/iy(k',f)) =P;^(k,t)(27r)35(k + k'). (22) 

The gravitational wave energy density power spectrum is 
then 


dpGW 

c?ln(fc) 


1 

327rG^ 


P^(k,f). 


(23) 


Strictly, Eq. (ITO)) applies only for detonations. For defla¬ 
grations, one should take into account the suppression of 
the tunnelling rate ahead of the bubble wall, where the 
fluid is heated by the release of latent heat. In this case 
we would expect P* ~ Cs//3. 

The important ratio /3(the transition rate relative 
to the Hubble rate) follows from simple considerations of 
the temperature of the transition . One can straight¬ 
forwardly argue that 


B. GW power spectrum from fluid and field 

The source of gravitational waves is the transverse 
traceless part of the spatial components of the energy- 
momentum tensor. Given that we will be removing 
the trace anyway, it suffices to consider a source tensor 
Tij = rfj -I- r/j, which is decomposed into fluid and field 
pieces according to 

4 = 4 = W2(e + p)ViV,. (24) 


^(tN) ^ 41n(mp/rN), (19) 

and that for tunnelling in a thermal effective potential 

0 


/? ^ 25'(tN) 

K ~ (i-Tn/Tc)' 


( 20 ) 


The physical metric perturbations are recovered in mo¬ 
mentum space by applying the projector onto transverse, 
traceless symmetric rank 2 tensors: 

Ajj.fei(k) = Pikik)Pji{k) - ipy(k)Pfcz(k) (25) 


with 


Hence, for a thermal electroweak-scale transition, the 
critical bubble action must be O(IO^), and the ratio /3/i7* 
must be at least O(IO^). 

A detailed non-perturbative evaluation of the bubble 
nucleation rate in the standard model electroweak the¬ 
ory is presented in Ref. [l^, using an unphysically small 
Higgs mass in order to ensure a first order phase transi¬ 
tion. In this case the critical bubble action was found to 
be ~ 90, and /3/i7* ~ 2 x lO'^. These are expected to be 
generic numbers for any first order thermal electroweak- 
scale transition. 


Pij{k) = S^j - kikj. (26) 

The particular solution for the gravitational wave is 
therefore 

%(k,t) = (167rG)Ay,fcj(k) J — —Tki{k,t'), 

(27) 

where we have assumed that the source vanishes for t' < 
0 . 

Using the fact that the fluid shear stress dominates the 
spatial parts of the energy-momentum tensor, we write 
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(28) 

Introducing the unequal time correlator (UETC) of the fluid shear stress iP [H, |4^ through 

A*j,fcz(k) (rf^ {k',t2)'^ = (fc, ti, t2)(27r)3(5(k + k') (29) 

and averaging over a period T, much longer than the periods of the gravitational waves of interest, we can write 

Pj^{k,t) = f —^n^(fc,ti,t2)- (30) 

Jo 2 


hk(t)hk'(^)) = (167rG)^y (itidt 2 cos[A:(t - ti)] cos[fc(t - t 2 )]Ay,fei (k) {^rf'’(k,ti)Tf"'(k',t 2 ) 


On dimensional grounds, we can write the UETC as 

n2(fc,U,t2):^ [(e + p)C7f]'T?n' (31) 

where e and p are the spatially averaged energy density 
and pressure; Uf is the root mean square fluid velocity, 
defined through 

(e + p)I7f = ^ y > (32) 

where V is the averaging volume; L{ is a characteristic 
length scale in the velocity field; and iP is a dimension¬ 
less function of k, ti and t 2 - In Ref. [s^ we estimated 

I 


that L{ would be the mean bubble separation, but we 
will not make that assumption yet. We will see that we 
can understand the numerical results better if we extract 
the scale directly from the fluid velocity field in the sim¬ 
ulations. 

We also assumed that the UETC would be a function 
of ti — O for times between the nucleation time In and 
the lifetime of the velocity perturbations Tv, and that 
there is no separate timescale in the function iP, apart 
from that generated from L[, the speed of sound Cs, and 
the speed of light. With these assumptions, we can write 
the dimensionless UETC as a function of kLf and z = 
k{ti — t 2 ), and the spectral density of h becomes 


Pf^{k,t) = 16'!TG{e+p)Uf 


tk 


■u?/ 


dz 


COs{z)^2 


W{kLi,z). 


(33) 


Note that one could follow through the same arguments for the scalar field, which would contribute in exactly an 
analogous manner 


pl{k,t) 


16 ;TG(e + p)C/^ 


tk 


-1 r3 


COS(Z) - 2, 


/ dz — 


(34) 


where 

{^ + P)ul = ^ J^d^xrfi, (35) 

is a characteristic scale in the scalar field configura¬ 
tion, and is the dimensionless unequal time correla¬ 
tor of the scalar field shear stress tensor. However, as 
explained in the introduction, the field contribution is 
negligible in most phase transitions. 

Hence, putting together (1^ and (1551) . we may write 
the gravitational wave energy density power spectrum as 


dpGW 

d\n[k) 


SttG 


-p)Uf 


tL{ 


{kLtf 

27r^ 


pQwikLf), (36) 


where 


PcwikLi) = ^J dz^^^ti\kU, z), (37) 


is a dimensionless spectral density for the gravitational 
waves. 

The gravitational wave power spectrum at time t can 
then be written 

PGW = (e+p)^Gf (<Lf)(87rGHGw), (38) 

where 

^GW = J *^^PGw(fcif) (39) 

is a dimensionless number. We see that the gravitational 
wave power spectrum grows linearly with time, for as 
long as the velocity perturbations are active, with a slope 
which depends on the square of the enthalpy density, the 
fourth power of the mean square fluid velocity, the fluid 
length scale, and a dimensionless number describing the 
fluid flow Hqw- 
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In principle, the value of flaw depends on the param¬ 
eters of the phase transition in dimensionless combina¬ 
tions, which we can expect to include the bubble wall 
speed Uw and the latent heat relative to the total en¬ 
ergy density ar- In Fig. 2 (bottom) of [sj, we plotted 


PGwl 


{e + p)'^UfLf 


against time. Noting that we have 


G = 1, the slope of the graph is SttIIgw- We found that 
Hgw was approximately constant, varying by no more 
than a factor 2, when we took the fluid s cale Li to be 
the mean bubble separation i?* = s/vJW^. Hence most 
of the dependence of the gravitational radiation energy 
density on the phase transition parameters is accounted 
for by the explicit factors in Eq. (l38l) . 


C. Integral scale 

The question of which scale to take for Lf affects the 
value of Hgwi and hence its variation between simula¬ 
tions. As mentioned above, in Ref. [s^ we took the scale 

I 


to be i?*, the average bubble separation. However, one 
could equally estimate the length scale from the veloc¬ 
ity field itself, and to this end we can use the following 
quantity (sometimes referred to as the integral scale) 


6 = 


1 

W) 


d^k 

(^ 


\k\-^Pv{k), 


(40) 


where is the RMS velocity. We will see that when 
the scale Lf is chosen to be the integral scale, the varia¬ 
tion in the parameter Hgw is reduced to about 10%. This 
emergence of Hgw as a quasi-universal constant for first 
order phase transitions with < 0.1 is an important 
result. 

One can also define an integral scale ^gw for the grav¬ 
itational wave energy density from its spectral density 
Pgw- We will also conhrm that the integral scale of the 
gravitational radiation is related to the integral scale of 
the velocity field, as one would expect. 


D. Dimensionless GW power spectrum parameter Ogw 


It is often useful to express the gravitational wave power spectrum as a fraction of the critical density, pc = SU^/SttG. 
Hence we are led to consider a dimensionless gravitational wave power spectrum 


dflGwjk, t) 
d\n{k) 


167rG(e -I- p)Uf 


2 tLf {kLif 
In 247r2 


PcwikLf), 


(41) 


where is the Hubble parameter at the time the bubbles are nucleated. Noting that the critical density is the 
energy density e, and denoting the lifetime of the source by Ty, we find that the dimensionless gravitational wave 
power spectrum during the radiation era can be written 


dHGw(fc) 

dln(fc) 


= 3{1 + wfUf{H^T^){H^Lf) 


{kLff 

2n 


PowikLf), 


(42) 


where w = p/e is the equation of state parameter. Inte¬ 
grating over wavenumber, we see that the total relative 
energy density is 

Hgw = 3(1 -k wfTf^{H,Tn(H^Lf)nGw- (43) 

E. Source lifetime 

It is clearly important for the calculation of the grav¬ 
itational wave energy density to calculate the lifetime of 
the source, the shear stress caused by the sound waves. 
We show in Appendix [A] that in an expanding universe, 
the shear stresses decay and decorrelate in such a way to 
make Ty precisely equal to the Hubble time. The shear 
stresses also decay due to the viscosity of the fluid at a 
scale-dependent rate. We should therefore estimate on 
which scales viscous damping time is smaller than Ty. 


For linear non-relativistic flows induced by sound 
waves (i.e. for velocity fields V/ which are purely lon¬ 
gitudinal) , viscosity adds a term of the form 

(^^s + Cb) (44) 

to the left hand side of Eq. (ini), where rjs is the shear 
viscosity and Cb is the bulk viscosity. For a plasma of 
relativistic particles in a gauge theory, the bulk viscosity 
is negligible compared to the shear viscosity , and the 
shear viscosity can be estimated as 

77s--r^/e'‘ln(l/e), (45) 

where e is the electromagnetic gauge coupling [5l| . Hence 
velocity perturbations of wavenumber k are damped as 
exp(—477sfc^t/3), and the lifetime due to viscous damping 
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of sound waves with wavelength R is 

T^(i?) ~ i?^e/r 7 s ln(l/e)i?^r. (46) 

Therefore, at a transition with temperature just below 
the critical temperature Tc, the viscous damping lifetime 
exceeds the Hubble time H~^ for all scales 



where we have neglected the logarithm of the gauge cou¬ 
pling. 

We will see in the next section that the scale of the 
fluid perturbations is set by the average separation of the 
nucleating bubbles i?*, and that the bubble separation at 
an electroweak-scale phase transition with any interesting 
degree of supercooling will satisfy this inequality. 

Hence for a first order transition at the electroweak 
scale - or even a few orders of magnitude above - the 
lifetime of the source of the gravitational waves is the 
Hubble time, 


r, = (48) 


F. Comparison to envelope approximation 


In the envelope approximation, the relative energy den¬ 
sity in gravitational waves is given by 


“gw — 


0.42 -k v\ 


2 2 2 
n a 


H, 


15 J (a -I-1)^ 


(49) 


where a is the ratio between the “vacuum” energy (de¬ 
fined below) and the radiation energy density in the sym¬ 
metric phase, K is the efficiency with which vacuum en¬ 
ergy is converted to kinetic energy, and /3 is the nucle- 
ation rate parameter also defined above. 

The vacuum energy Vq is defined in Ref. [2^ from the 
trace anomaly. 


e = e-3p, (50) 

as a quarter of the difference between the symmetric and 
broken phases: 


Vo = \{0s-Oh)- (51) 

In our convention, the trace anomaly vanishes in the sym¬ 
metric phase, and in the broken phase is 

0b = -T^H(<^b,T) + 4H(0b,T), (52) 

where is the value of (f) in equilibrium in the broken 
phase at temperature T. In the conventions of [1^, the 
trace anomaly vanishes in the broken phase, and in the 


symmetric phase is equal and opposite to (1521) . Hence for 
our thermal potential m the parameter a is 


a = = ^-7 (-T—V((l)h,T) - I/(())b,T)^ . 

3aT4 3ar4 ^4 dT ’ 7 

(53) 

The efficiency parameter is defined from the average fluid 
kinetic energy density (l32|) as 

Therefore 

(l+u;)I7f = -^. (55) 

1 -I- a 

The factor of 1 -I- a in the denominator of the right hand 
side comes from the fact that we are dividing by the av¬ 
erage total energy density in the symmetric phase, which 
is 3aT^ + Vb in the conventions of [2^ . 

Note that na is conventionally estimated analytically 
from the radial fluid velocity around an isolated expand¬ 
ing bubble u(r, t), where r is the distance from the centre 
of the bubble, and t is the time since nucleation [2^, [2^. 
At large times, the radial fluid velocity is a function of a 
scaling variable ^ = r/t, rather than r and t separately. 
The ratio of the kinetic energy density to the total energy 
density can then be estimated as 

[ d^f{e + p)W^v^{0- (56) 

vie J 

We will compare this estimate to the numerically ob- 
_2 

tained {1 + 'w)Uf in the results section, finding good 
agreement. 

In order to compare our expression for the gravita¬ 
tional wave energy density (1431) with the envelope approx¬ 
imation formula (HHl) . we estimate the fluid flow scale L{ 
as the bubble separation scale i?*, which in turn is related 
to the nucleation rate parameter by Eq. dig. 

Hence the ratio between the gravitational wave energy 
density generated acoustically and in the envelope ap¬ 
proximation if0 


Ogw 


ft 


ea 

GW 


3(87r)3HGW ,n ^ 
0.1Iu2(0.42 + u2) 


(57) 


Given that the ratio (lATl) is smallest for tiw = 1, and that 
the lifetime of the sound waves is approximately 
(see Section fill El) , we can estimate that 


Hgw 


> 60Hgw 


Qea VJVV ^ • 


(58) 


^ Note that for a deflagration, if tunnelling is suppressed behind 
the shock wave, the ratio is boosted by a factor Cs/uw - see 
the discussion after Eq. (II 8 II . 















We will see from our numerical simulations that ^gw ~ 
0.04. The ratio jS/H^ was discussed in Section HlBl and 
shown to be at least O(IO^), and possibly significantly 
greater if there is only small supercooling. We conclude 
that the energy density in acoustically generated gravi¬ 
tational waves is at least two orders of magnitude greater 
than the envelope approximation suggests. 


IV. NUMERICAL SIMULATIONS 

A. Methods 

Our numerical methods are a development of those first 
used in this context to study the case of isolated bubbles 
in Ref. [s^. In that paper a spherically symmetric bub¬ 
ble was assumed. Here we extend those simulations to 
a full 3-|-1-dimensional simulation volume. In addition, 
we couple the linearised stress-energy tensor to perturba¬ 
tions around a flat metric, to measure the gravitational 
wave power produced by the simulation. 


1. Coupled field-fluid system 

The coupled hydro-scalar equations, outlined above, 
can be treated quite easily using standard numerical tech¬ 
niques. The scalar field is evolved with the leapfrog (Ver- 
let) algorithm, while standard operator splitting methods 
are used for the fluid . These are equivalent to numer¬ 
ically integrating the equations of motion given above. 

Although the full details of how to implement relativis¬ 
tic hydrodynamics is beyond the scope of this paper, it 
is instructive to consider how the quantities are laid out 
on the lattice both in the spatial and temporal directions 
(see Fig. [T]). Furthermore, for good energy conservation 
it is essential that the discretised version of the damping 
term couple the field and fluid quantities at equal times 
during the simulation. 

We have tested the results of our simulations against 
changing timestep (as well as the lattice spacing); see the 
following section. 

As our simulations do not run for sufficiently long to 
develop strong shocks (indeed, we choose our lattice spac¬ 
ing parameters such that the fluid velocity profile is al¬ 
ways resolved by several Sx), the simulations presented 
in this paper do not involve any artificial viscosity. The 
importance of an artificial viscosity term was previously 
studied using l-|-I-dimensional simulations of two collid¬ 
ing bubbles. 


2. Metric perturbations 

Our principal observables are the energy density and 
power spectrum of the gravitational waves. The goal 
of our simulations is to compute the power per unit 



e, T, p, 

(j), TT 


u, z 

y 

z] 

* -> X 

Vy 



e, T, p, (j) 


U, V,Z 

t 

TT 



U 


FIG. 1. Layout of quantities simulated. The positions of 
quantities related to simulating an ideal relativistic fluid are 
standard [^. Because the field and fluid are coupled to¬ 
gether, it is important that the scalar field (j> and its conjugate 
momentum tt are correctly centred. We take 0 to reside in 
zones (like pressure, temperature, etc.), so that no centring is 
required to compute, for example, the equation of state. 

logarithmic frequency interval in gravitational waves 
dpGw(^)/dlnA:, and the total energy density PGW- 

Perturbations of the metric are sourced by transverse- 
traceless part of the stress-energy tensor Hy 

hy - = IGTrCHy . (59) 

Obtaining Hy from involves a projection in momen¬ 
tum space. Therefore, evolving hy (whether in momen¬ 
tum space or position space) would involve Fourier trans¬ 
forms at each timestep. As we go to large volumes, the 
execution time of fast Fourier transforms (FFTs) scales 
as 0{NlogN), while there are few optimised FFT codes 
that offer domain decomposition in more than one dimen¬ 
sion. It is therefore vital that the number of steps requir¬ 
ing Fourier transforms be minimised, to yield a scalable 
simulation. 

Our approach is to evolve the unprojected equation of 
motion in real space 

iiij - = WirGirf^ + 4 ) 


(60) 
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(384/r^J^ volume, at time 1000/T^ 



FIG. 2. Single bubble test simulation, with the correlation 
length I shown as an indication of the wall width. Only the 
fluid source is shown here; discretisation errors for the field 
source are the same or smaller. There is good agreement 
between 1/L and l/l, as desired. 

where Uij is an auxiliary tensor and the sources are de¬ 
fined in Eq. (1241) . Only when we wish to recover the 
metric perturbations hij do we Fourier transform Uij and 
project out the transverse-traceless components through 

(k) — (kjrt/TTi(t, k), ( 61 ) 

where the projector is defined in Eq. (|25|) . We evolve 
Eq. (1601) using a leapfrog algorithm in a similar manner 
to the scalar field. 

Note that we choose the units of the code such that 
the critical temperature Tc = 1 and the gravitational 
constant G = 1. 


B. Tests 

Our basic tests principally involve varying the lattice 
spacing and timestep independently, on simulations of 
a single bubble colliding with itself in a small periodic 
box. These allow us to test that the simulations perform 
accurately between length scales l/R* and 1/i. Longer 
distances do not need to be tested, and in any case i?* is 
set by the box size L in these tests. 

1. Changing the lattice spacing 

We performed tests on the effect of changing the lat¬ 
tice spacing using the self-collision of a single bubble 
in a cubic box, in relatively modest volumes (with pa¬ 
rameters given in the following section). We considered 
Sx = 0.5/Tc, Sx = 1/Tc, Sx = 2/Tc and Sx = 4/Tc. 
We notice no significant difference between these choices 
until Sx = 4/Tc. 

It is worth mentioning that, even for an isolated bubble 
which would (in continuum) have vanishing quadrupole 
moment and hence not source gravitational waves, the 


lattice discretisation breaks the spherical symmetry and 
results in a small amount of gravitational wave produc¬ 
tion. This power goes to zero as (Ja;)'^ for both the field 
and the fluid sources. After collision, however, agreement 
is very good with relative differences of at most 7% for 
k < 1/i between Sx = 1/Tc and Sx = 2/Tc for the fluid 
source. Furthermore, at higher momenta there are only 
0(1) differences between these two choices, consistent 
across seven orders of magnitude. This is surprisingly 
good given the relatively coarse wall width and the com¬ 
plicated microphysics. Similarly, discrepancies between 
Sx = 0.5/Tc and Sx = 1/Tc at late times were at worst 
2% for k < l/i\ see Figure [2j Discretisation errors were 
always less severe for the field source than for the fluid 
source. 

In summary, we note no significant sensitivity to lattice 
spacing so long as it is kept well below the scalar field wall 
width. 

While our previous work used simulations with Sx = 
1/Tc, we use a lattice spacing of Sx = 2/Tc in the 
present paper. The inferred discrepancies are demonstra¬ 
bly smaller than 10%, and the doubling of the accessible 
dynamic range that this allows is very useful. 

2. Changing the timestep 

With Sx = 2/Tc having been chosen, we varied the 
timestep to explore the effect of inaccuracies in our evo¬ 
lution algorithm. There is agreement at the 1% level or 
better for all k < 1/i and 5% or better up to /c ~ 0.5 
(all points plotted on Fig. [2]) as we varied St between 
Q.2/Tc, 0.1/Tc and 0.05/Tc, in the same single-bubble 
tests for the fluid power spectrum described above. We 
use St = 0.1 for the remainder of the paper, although 
we could probably have achieved acceptable results with 
St = 0.2. 

In the present paper our simulation durations are typ¬ 
ically the same order of magnitude as one light-crossing 
time, and rather less than one sound-crossing time. This 
means, in particular, that the production of gravitational 
radiation by acoustic waves (or by scalar radiation, which 
is in any case heavily damped) is not likely to be affected 
by signals propagating around-the-lattice. 

C. Parameter choices 

We use the same parameters for the potential as in our 
previous paper. The exact values of these parameters are 
not particularly important: it is the latent heat and the 
wall velocity which mainly determine the gravitational 
wave power spectrum. Our aim in the present paper is 
to develop the ideas underlying our previous paper as 
well as the spherical studies carried out earlier, and so 
we work with the same parameters as before. 

No attempt is made in the present paper to look at 
strong fluid flows or fast ‘runaway’ bubble walls. We 
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Weak 

Weak (scaled) 

Intermediate 

To/Tc 

1/V2 

1/V2 

1/V2 

7 

1/18 

4/18 

2/18 

A 

^/72 

^/9 

^/72 

A 

10/648 

160/648 

5/648 


9/40 

9/40 

9/5 


1/10 

1/20 

4^2/10 

ITc 

6 

3 

6/^2 

Ttiin 

0.86 

0.86 

0.8 


0.010 

0.010 

0.084 

a 

0.0046 

0.0046 

0.050 

RcTc 

16 

8.1 

8.6 


TABLE I. Scalar potential parameters ©, nucleation tem¬ 
perature Tn, phase transition parameters (na , transition 
strength parameters and (isa, and critical bubble radii 
IMJ for our simulations. 


leave these harder topics for future work and instead 
seek to comprehensively explain the generation of gravi¬ 
tational waves by more gentle phase transitions (aTj., ^ 
0 . 1 ). 

We discuss bubble nucleation further in the next sec¬ 
tion but note that we nucleate all of our bubbles simul¬ 
taneously in the present work. 

Given 6x = 2/Tc and a simulation size of 2400^ points, 
our physical simulation volume is (4800/Tc)^ for all the 
results presented in this paper. 


D. Initial conditions 


At the start of our simulation, we nucleate a con¬ 
trollable number of bubbles, which was usually 0(1000) 
(yielding bubbles of average collision radius slightly 
larger than in Ref. M) , but was as small as 37 or in 
one case as large as 32558. These have a Gaussian scalar 
field prohle. This profile is initially at rest, meaning that 
the conjugate momentum, and also the fluid velocity are 
zero in the vicinity of the bubble. We ensure that all 
the initially nucleated bubbles are well separated at the 
start of the simulation. For runs with the same num¬ 
ber of bubbles but different wall velocities, all bubbles 
are nucleated at the same positions, but from testing we 
found that even 37 bubbles was enough to remove any 
discernible dependence on the initial bubble positions. 

The critical bubble radius can be computed from the 
surface tension a and the difference in potential energy at 
Tn from the thin-wall formula (noting that the potential 
energy in the symmetric phase is zero) 


2cr 

-R((/)b,rN)' 


(62) 


Values of Rc for our simulations are shown in Table ID 
Rather than find the critical bubble profile exactly, we 


use a spherically symmetric Gaussian field profile 

= (f>b exp(-r^/2Rc)- (63) 

This is rather broad, and therefore sufficiently large com¬ 
pared to the true critical bubble profile to ensure that the 
bubbles reliably expand despite lattice effects. 

The bubbles are sufficiently large that they immedi¬ 
ately start growing, driven by the pressure difference 
between the interior and the exterior. The scalar field 
quickly settles into a kink-like configuration, interpolat¬ 
ing between the metastable and stable minima over a dis¬ 
tance of order £, the correlation length of the scalar field 
(see Table H] for the values the correlation length takes). 
For the scalar field dynamics to be valid we must have a 
lattice spacing that resolves the wall width (see previous 
section), which places an upper limit on the physical sim¬ 
ulation volume possible for a given amount of computer 
memory. 

In this paper, the bubbles are nucleated simultane¬ 
ously. Nucleating at a single time helps to ensure clear 
scale separation in the limited dynamic range avail¬ 
able to our numerical simulations, although it does pro¬ 
duce oscillatory patterns in the resulting power spectrum 
(we cover the case of unequal nucleation times in Ap¬ 
pendix [B]). We could in principle recover the power spec¬ 
trum produced by bubbles of all different sizes by a linear 
superposition of the resulting power spectra, weighted by 
the bubble size distribution. 

Once nucleated, the bubbles grow, and the fluid ap¬ 
proaches a characteristic radial velocity distribution, 
which is a function of ^ = r/t, where r is the distance 
from the centre of the bubble, and t is the time since nu¬ 
cleation (see Fig. [S]). The form of this function depends 
on the bubble wall velocity , and we will refer to it as 
the scaling profile. The rate of approach to the fluid scal¬ 
ing profile is generally much slower than the relaxation 
of the scalar field. In a true electroweak phase transition 
the bubble size at collision is many orders of magnitude 
larger than the bubble size at nucleation, giving a lot 
of time for the radial velocity distribution to reach its 
asymptotic profile. 

In our numerical simulations, the ratio of the bubble 
size at collision, to the bubble size at nucleation (rs 
Rc) is at most 90 for W = 37, i = 16 and as small as 9.4 
for W = 32558, £ = 16 (our simulation parameters are 
outlined in Tables U and |lTl) . We should therefore be alert 
to the fact that the fluid has definitely not settled down 
to its final scaling profile in a collision. One can test 
for the effect of a non-scaling fluid profile by repeating 
simulations with fewer bubbles, so that there is a longer 
time before collision. We have carried out simulations 
such that R varies by around a factor of three in the two 
sets of simulations for which we present plots, and by as 
much as a factor of ten in our full set of simulations for 
this paper. 
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FIG. 3. Comparison of radial fluid velocity profiles for sim¬ 
ulations at approximate collision times for = 1000 {t = 
500/Tc', gray), A^b = 37 (t = lOOO/Tc; black) and at late times 
(dashed red). 


Type 

V/E 

«w 

-^b U max 

TT^ 

^ f, max 

Cf, end^c 

00 

0 

Weak 

0.06 

0.83 

988 

0.0052 

0.00037 

351 

0.88 




125 

0.0052 

0.00028 

649 

0.84 


0.1 

0.68 

988 

0.0084 

0.00036 

244 

0.73 




125 

0.0082 

0.00026 

451 

0.71 




37 

0.0080 

0.00021 

644 

0.60 


0.121 

0.59 

988 

0.0116 

0.00052 

182 

0.69 


0.15 

0.54 

988 

0.0102 

0.00037 

230 

0.54 




37 

0.0120 

0.00025 

428 

0.80 


0.2 

0.44 

32558 

0.0059 

0.00047 

136 

0.97 




988 

0.0073 

0.00031 

368 

0.70 




125 

0.0075 

0.00023 

613 

0.86 




37 

0.0078 

0.00019 

942 

0.70 


0.4 

0.24 

988 

0.0036 

0.00049 

756 

0.86 

Wk. (sc.) 

0.4 

0.44 

988 

0.0075 

0.00029 

365 

0.81 

Interm. 

0.4 

0.44 

988 

0.0595 

0.00328 

485 

1.04 


TABLE II. Simulation parameters rj (field-fluid coupling), Nh 
(number of bubbles nucleated), with the resulting bubble wall 
speed Uw, the maximum fluid RMS velocity Ut^ max, the max¬ 
imum contribution of transverse fluid motion [/f maxi Ihe in¬ 
tegral scale of the fluid end, and the scaled slope parameter 
for the growth of the gravitational wave energy density flew. 
The potential parameters and derived quantities for each type 
- “weak”, “weak scaled” and “intermediate” - are given in 
Table H 


E. Scaling 

A cosmological first order phase transition is a multi¬ 
scale problem, with length scales varying from the mi¬ 
croscopic (1/T, bubble wall thickness) up to the Hubble 
scale, a range spanning 17 orders of magnitude at the 
electroweak scale. The typical bubble sizes at collision 
time are somewhere between these scales, depending on 
the metastability of the high temperature phase. It is of 
course impossible to include all of these scales in a sin¬ 


gle numerical simulation, where scale hierarchies only of 
order 10^ are achievable. To obtain a stable numerical 
description of the bubble wall, the wall thickness has to 
span a few lattice units (denoted by 6x). In order to have 
collisions within the simulation volume, this restricts the 
bubble separation to be unphysically small. 

However, this restriction can be relaxed, at least 
partly: we expect that the dynamics of the bubble 
growth, collisions and the subsequent generation of gravi¬ 
tational waves are mostly determined by the “bulk” ther¬ 
modynamics (e, p, latent heat C) and the friction param¬ 
eter r]^ but not by microscopic details of the bubble wall 
(surface tension a, wall thickness i). Dimensionally, it is 
clear that the contribution from quantities proportional 
to bubble volume (e.g. latent heat) will dominate over 
quantities proportional to the area of the bubbles when 
the bubble radius is large enough. 

This motivates us to search for a way to modify the 
equations of motion (nnn-dni) so that we could simulate 
bubbles which are significantly larger than the micro¬ 
scopic length scale, while preserving the bulk thermo¬ 
dynamics of bubble expansion while possibly sacrificing 
the properties of the bubble wall. Indeed, this can be 
achieved with the following simple rescaling of the pa¬ 
rameters and fields: 

7 —>■ r^7, A —>■ A —>■ r^A, rj —>■ ryy, 

—>■ r~^(j){rx), V^{x) —>■ lA*(rx), (64) 

E{x) —)► E{rx), Zi{x) —)> Zi{rx), 

Here r is a dimensionless scaling factor, and x = (x, t). 
Clearly, the equations of motion (nnii-dnD remain valid. 
The crucial feature of the scaling is that the potential 
remains constant, V{(j),T) = I4caied(?'~^(/', T), indicating 
that the bulk quantities Tc, C and also e(T, </>) and p(T, (/>) 
remain invariant, as desired. However, the surface ten¬ 
sion and wall thickness scale as cr —>■ r“^(T and £ —>■ 

In effect the scaling stretches the field configuration by a 
factor of r~^ in spatial and temporal directions. 

We note that in spite of the non-trivial scaling of the 
friction parameter 77 , the total frictional force imparted 
on the moving bubble wall does not change: it is ob¬ 
tained by integrating the 77 -terms in Eqs. (nnD-dnD over 
the bubble wall thickness, which is scaled by a factor of 
r~^. 

What does the rescaling gain us? It is straightforward 
to see that the lattice implementation of the equations 
of motion (fTUl) - (fT^ do not change (in lattice units) un¬ 
der scaling (IMl) . provided that the lattice spacing is also 
scaled as Sx —>■ r~^ 6x. This implies that a single lat¬ 
tice simulation exactly corresponds to a whole family of 
results, given by the scaling with r. All of them have 
the same bulk thermodynamical properties. Thus, pro¬ 
vided that the detailed bubble wall properties are not 
important for bubble collisions and gravitational wave 
generation, we can take a simulation run where bubbles 
have been nucleated at specific locations, and rescale it 
to the desired physical bubble separation scale. 
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We can test the assumption that the surface properties 
are not important by comparing results from simulations 
which differ only in surface tension and wall thickness. 
This can be achieved by applying the scaling (IMl) to the 
parameters of the theory, but leaving the lattice spacing 
constant. A set of parameters for unsealed (r = 1) and 
scaled (r = 2) runs are shown in TableHl The surface ten¬ 
sion and the bubble wall thickness have been halved in 
the scaled simulation. The results of the test are shown 
in Table m here the scaled run is done with rj/Tc = 0.4 
and 988 bubbles, which can be compared with unsealed 
t]/Tc = 0.2, 988 bubble results. In both simulations the 
bubbles are nucleated at identical times and locations. 
The numerical results match well within uncertainties of 
the measurements, supporting our assumption that the 
surface properties of the scalar field profile are unimpor¬ 
tant. 


V. NUMERICAL RESULTS 

In this section we present the main results from a 
campaign of numerical simulations, whose parameters 
are given in Table |TT1 As mentioned before, our simu¬ 
lations were carried out in a volume (4800/Tc)^. A set 
of representative slices through the simulation are shown 
in Fig. IT] Our main results are derived from the set 
of simulations with latent heat to thermal energy ratio 
q;tn — 0 . 01 , which we characterise as a “weak” tran¬ 
sition. Values of the friction parameter ry were chosen 
to give bubble growth proceeding by both detonation 
and deflagration, as well as one simulation tuned to the 
Jouguet case, where the bubble wall moves at the speed 
of sound. We have one “intermediate” strength transi¬ 
tion, with q;tn — 0.1 where 77 is chosen to give the same 
wall speed as the weak transition with ry/Tc = 0.2. The 
“weak scaled” transition (employing the scaling of the 
previous section) is discussed later. 

When plotting graphs, we focus on three representa¬ 
tive cases, where the field-fluid coupling is rj/Tc = 0 . 1 , 
rj/Tc = 0.15, and rj/Tc = 0.2, and the bubble wall speed 
is supersonic (v^ = 0.83), just subsonic (t^w = 0.54), and 
subsonic (fw = 0.44). A complete set of graphs can be 
found in the supplementary material [53 |. 

Our understanding of the transition developed in Sec¬ 
tion [ITT] shows that the important quantities for the over¬ 
all gravitational wave energy density are the RMS fluid 
velocity Uf and the fluid velocity scale if, and that the 
gravitational wave power spectrum is only indirectly de¬ 
pendent on the strength of the transition and the pa¬ 
rameters in the potential. Indeed, the gravitational wave 
power spectrum should be the same for parameters which 
give the same q;tn — 0.01 (keeping the wall velocity con¬ 
stant). We can use the “weak scaled” run of Table HIl to 
test this statement, where otn is constant but the scalar 
bubble wall thickness is halved. We test the effect of the 
strength of the transition with the “intermediate” run of 
Table im 


We track the progress of the transition through the 
time evolution of the two quantities C/^ and Uf defined 
in Eqs. (1321) . (1351) . We recall that the squares of these 
quantities give an estimate of the size of the shear stresses 
of the field and the fluid relative to the background fluid 
enthalpy density, and that the Uf tends to the RMS fluid 
velocity for l/f <C 1. We also note that the fraction of the 
fluid velocity power coming from rotational modes, U f , 
is very small leading us to conclude that rotational fluid 
modes are not important in this system; we discuss this 
in more detail in Appendix [Cl 

We see from Fig. [5] that U^ grows and decays with 
the total surface area of the bubbles of the new phase, 
while the mean fluid velocity grows with the volume of 
the bubbles, and then stays constant once the bubbles 
have mergecd. This allows us to identify distinct phases 
of the transition; the collision phase, where U ^ grows 
and decays; and the subsequent acoustic phase where Uf 
is approximately constant, and U^ vanishes. 

A. Length scales 

The analysis of Section Hill shows that the length scale 
of the velocity flow is an important determinant of the 
gravitational wave power spectrum. In Fig. [ 6 ] we show 
the integral scales for the velocity and the gravitational 
radiation for runs at = 37 and TVb = 988. During 
the collision phase, the bubbles expand and overlap, and 
hence the scales of the velocity field and the resulting 
gravitational radiation grow linearly in time. The gravi¬ 
tational radiation length scale is 2-3 times that of the ve¬ 
locity field. The scale of the velocity field stops growing 
as the bubbles collide and the scalar field decays to the 
vacuum, and stays constant during the acoustic phase. 
The scale imprinted on the gravitational radiation dur¬ 
ing the acoustic phase is close to that of the velocity field. 

B. Velocity profile 

Given the discussions on initial bubble sizes in Sec¬ 
tion [IV 0 it is important to bear in mind that the bub¬ 
bles in our simulations expand in size by a factor of only 
around 10 - 100 , which is many orders of magnitude less 
than in a real phase transition. One practical effect is 
that the profile of the velocity field around the bubbles 
does not reach its asymptotic scaling form, which can 
be expressed in terms of the previously introduced ra¬ 
tio f = r/t. In Fig. Owe showed the velocity profiles 
for the weak transition at rj/Tc = 0.1, rj/Tc = 0.15, and 
rj/Tc = 0.2, after times t = 500/Tc and t = lOOO/Tc- 


® We have no explicit viscosity, and the slight decreasing trend in 
some measurements of f/f arises from the well-known numerical 
viscosity of donor-cell advection, i^num — UfSx. 
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FIG. 4. Slices of fluid kinetic energy density EjT^ at t = 500Tc t = 1000To ^ and t = 1500Tc ^ respectively, for the 
rj/Tc = 0.15, = 988 simulation. 






FIG. 5. Root mean square fluid velocity U{ and root mean 
square scalar gradients for A^b = 988 (top row) and A^'b = 
37 (bottom row). 






FIG. 6. Plot of integral scales and ^gw associated with the 
fluid and gravitational wave power spectrum for simulations 
of interest for Nh = 988 (top) and A'^b = 37 (bottom). Note 
the different y-axis plotting scales. 


These are approximately when most bubble collisions are 
happening, in the A^b = 1000 and = 37 runs respec¬ 
tively. 

We see that, at collision, the velocity profiles are qual¬ 
itatively similar to their asymptotic forms in amplitude 
and shape, but differ in detail. In particular, the peak 
velocities are lower. This is particularly noticeable at 
the earlier time. We would therefore expect the RMS 


velocities Uf measured in the simulations to be under¬ 
estimates. As the gravitational wave power spectrum 
depends on the fourth power of U{, this is a significant 
source of uncertainty in deriving accurate predictions for 
the gravitational wave power spectrum. 

These considerations are tested in Table IIIII where 
we compare our mean square fluid velocity parameter 

Uf with which should be equal according to the 
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V/Tc 

Ww U{ 

\/ jKldOl 

4 

0.06 

0.83 0.0052 

0.0056 

0.0063 

0.1 

0.68 0.0084 

0.0085 

0.0121 

0.121 

0.59 0.0116 

0.0146 

0.0192 

0.15 

0.54 0.0102 

0.0103 

0.0100 

0.2 

0.44 0.0073 

0.0066 

0.0065 

0.4 

0.24 0.0036 

0.0033 

0.0036 


TABLE III. Simulation parameters rj (field-fluid coupling), 
with the resulting bubble wall speed Uw, fluid RMS velocity 
U{, for weak transitions with Nh = 988, and the equivalent 
quantity yIko/S appearing in the envelope approximation 
(see Eq. [SS]). The efficiency parameter k is estimated in two 
ways: Kiaa is estimated from the numerical spherically sym¬ 
metric ID fluid profiles at t = 1000/Tc, while kbsp comes from 
the function k(uw, a) given in the Appendix of Ref. [^ . using 
Uw extracted from spherical ID simulations at f = 1000/Tc. 

discussion in Section IIIIFI In the table, the efficiency 
parameter n is estimated in two ways: Kidcx is esti¬ 
mated from integrating the numerical ID fluid profiles at 
t = 1000/Tc, while kesp comes from the function k(uwj ct) 
given in the Appendix of Ref. [2^, using Uw extracted 
from ID simulations at t = 1000/Tc. As can be seen, Uf 
from the 3D simulations compares reasonably well to its 
estimate extracted from the ID numerical profiles around 
the time of bubble collision, while the theoretical values 
are somewhat higher. It is remarkable that such a sim¬ 
ple model for the mean square velocity, which omits all 
details of the bubble collisions, does so well. 

C. Power spectra 

In Figs. [7] and E] we show velocity and gravitational 
wave power spectra at various times through the simula¬ 
tions, for weak transitions with rj/Tc = 0.1, rj/Tc = 0.15, 
and rj/Tc = 0.2, where the bubble wall speed is super¬ 
sonic (uw = 0.83), just subsonic (uw = 0.54), and sub¬ 
sonic (uw = 0.44). The same potential and fluid-field 
parameters are run with A^b = 988 and A^b = 37 bubbles, 
to show the effect of allowing a greater time for the fluid 
velocity around the expanding bubbles to approach their 
scaling profiles. The power spectra develop in character¬ 
istic ways in the different phases of the transition, and 
one can see that if the simulation is stopped too early, 
a misleading impression of the power spectrum will be 
obtained. 


1. Collision phase 

Looking first at the velocity power spectra, the most 
striking feature is their periodic modulation. This is not 
a physical feature, and is due to the bubbles being nucle¬ 
ated all at the same time. We have checked that spread¬ 
ing the nucleation times reduces this modulation, and it 


is not expected to be a feature of the velocity power spec¬ 
trum of a realistic bubble nucleation distribution in the 
infinite volume limit. In Appendix [B] we show the effect 
of allowing nucleation over a time of about 200/Tc. 

Once the fluid shells of the nearest pair of bubbles be¬ 
gin to overlap, gravitational waves are generated, at a 
scale controlled by the size of the bubbles. The overlap 
of the fluid shells is quickly followed by the collisions of 
the bubble walls, and gravitational radiation is generated 
by the scalar field as well. The bubbles continue to grow 
and to collide, and as a result the length scales of the 
velocity field and the gravitational radiation get larger 
(see Fig. | 6 ]). This effect can be seen in the power spec¬ 
tra, where the curves show a peak moving up and to the 
left with time. 

In our simulations there is generally more energy in the 
scalar field than in the fluid to begin with, and so the 
gravitational radiation from the scalar field dominates 
the early phases (see Fig. 0 . However, when scaled to 
a real deflagration or detonation in the early universe, 
most of the latent heat of the transition goes to the fluid, 
and the radiation from the scalar field can be neglected. 
It is only in the case of a runaway bubble wall that the 
scalar field takes most of the latent heat. We discuss the 
scaling to real transitions in Section IVDl and we plan to 
study runaway transitions elsewhere. 

However, it is interesting to study the difference be¬ 
tween a fluid-only gravitational wave power spectrum, 
and one sourced by both fluid and field (Fig. 0. There 
one sees evidence of a power spectrum in arising from 
the scalar field during the collision phase (solid lines), 
which is later dominated by the gravitational waves from 
the fluid. As the scalar field energy density is confined 
to a thin shell, it is reasonable to suppose that its con¬ 
tribution can be adequately computed in the envelope 
approximation in the collision phase. We will investigate 
this conjecture elsewhere. 


2. Acoustic phase 

Eventually, the low-temperature phase spreads 
throughout the volume, the scalar field domain walls 
disappear, and fluid velocity perturbations are left 
behind. We call this the acoustic phase of the transition, 
as the fluid perturbations are primarily compressive (lon¬ 
gitudinal) modes (see Appendix [0. During the acoustic 
phase, the length scale of the fluid perturbations and 
the gravitational waves remains constant. 

For the simulations with fewer bubbles (Aj^ = 37), we 
see from the lower row of Fig. [7] that the envelope of 
the velocity power spectrum has an approximate power- 
law envelope beyond the peak. This power-law envelope 
is also visible at A^b = 988 at r]/Tc = 0.2, where the 
bubble wall speed is lower. In both cases at 77 /Tc = 0.2, 
the bubbles expand longer before collision, and we expect 
the velocity field to be closer to the asymptotic form. We 
note that the power law is approximately k~^ at rj/Tc = 
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FIG. 7. Velocity power spectra, for weak transitions, at rj/Tc = 0.1, 0.15 and 0.2 (vw = 0.83, 0.54 and 0.44) for Vb = 988 (top 
row) and Nh = 37 (bottom row). The large oscillations are due to all the bubbles being nucleated at exactly the same time. As 
in Fig. |8] we note that the scales are standardised for all the plots, but that the phase transition has not necessarily finished 
by 2500/Tc, the time of the latest curve. 


0.2, and appears steeper for lower couplings. However, 
we are not confident that we have reached the asymptotic 
form for bubble wall speeds above = 0.44. 

At low wavenumbers, the velocity power spectrum be¬ 
haves as a power of fc, and arguments based on the analyt- 
icity properties of the Fourier transform of a longitudinal 
vector field in Ref. show that it should go as fc®. This 
is just visible in the first few bins of the simulations with 
fVb = 988. Larger simulations are required to properly 
check the long-wavelength behaviour. 

We see that the gravitational wave power spectrum 
grows linearly with time in the acoustic phase, main¬ 
taining its shape, except at the lowest wavenumbers. A 
power-law behaviour can be seen emerging beyond the 
peak, especially in the simulations with = 37. The 
power-law is approximately k~^ for the weak deflagra¬ 
tions at rj/Tc = 0.2 and rj/Tc = 0.4 (the power spectra 
for the latter can be found in the supplementary material 
for this work) for both = 988 and W = 37, which 
gives us confidence that we are close to the true power 
law. However, a power-law can be seen only for W = 37 
for Tj < 0.15. Without further simulations at larger i?* 
we cannot properly determine the long-wavelength be¬ 
haviour of the growing acoustic phase power spectrum in 
these cases. 

We argued earlier that the gravitational wave den¬ 
sity parameter Hqw = PGw/e is proportional to L[, the 
fluid velocity length scale, and the square of the volume- 
averaged fluid energy density {e+p)‘^Uf. We plot the 


scaled gravitational wave energy density in Figure 1101 
This plot shows nicely parallel, linear growth of gravi¬ 
tational wave power when rescaled by these quantities 
at late times. The coincidence of the slopes is greatly 
improved over the equivalent figure in Ref. [s^, thanks 
to the larger simulation volumes, longer run times, and 
above all the replacement of the average bubble separa¬ 
tion at collision by the fluid integral scale. 

The improved coincidence of the slopes is one of the 
major results of the paper. It establishes the existence 
of an 0(1) parameter 87rHGw for a wide range of rele¬ 
vant transitions, and shows that the gravitational wave 
energy density from a phase transition can be understood 
in terms of simple features of the velocity field created 
by the dynamics of the bubble collision. 


D. Extrapolating to a real phase transition 

Our simulations are necessarily limited in volume, du¬ 
ration, and resolution. We now discuss how they can be 
extrapolated to the real universe. In particular we would 
like to extrapolate the gravitational radiation power spec¬ 
trum, expressed as a fraction of the critical density. 

There are three physical length scales in the system: 
the average bubble separation i?*, the size of the ini¬ 
tial bubble of the broken phase Rc, and the bubble wall 
width £. They are all set by the dimensional scale of the 
effective potential, which one can chose to be the critical 
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Ti=0.1 r,A'|_=988 



T)=0.1 r,A'|,=37 



11=0.15 r,iV^=988 



r|=0.15 r,iV|_=37 
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FIG. 8. Gravitational wave power spectra, for weak transitions, at rj/Tc = 0.1, 0.15 and 0.2 (vw = 0.83, 0.54 and 0.44) for 
A^b = 988 (top row) and A^b = 37 (bottom row). Note that the axes and time intervals are the same for all plots, which means 
that in some cases the latest (2500/Tc) curve is from before the completion of the phase transition. 



FIG. 9. Power spectra for r^/Tc = 0.2, comparing flnid-only 
(dashed) and total (solid) GW power at intervals of bOO/Tc. 
The power laws visible in the ‘total GW power’ case are dom¬ 
inated by the gradient energy of the scalar field. This sonrce, 
however, is short-lived. We conjectnre that it can be calcn- 
lated by means of the envelope approximation. 

temperature Tc, and various combinations of the dimen¬ 
sionless couplings 7 , A and A. In a real transition, the 
average bubble separation is much larger than the wall 
width because of the exponential factor in the tunnelling 
rate, whose argument is set by the ratio of the energy of 
the critical bubble (see Ref. [^) to the critical temper¬ 
ature. This is generally a large number. 

There are also two physical time scales to consider: 
the lifetime of the fluid flow Tv, and the duration of the 
phase transition, which is of order /3“^, the inverse of 
the tunnelling rate parameter uni). The duration of the 



FIG. 10. Time series of pG-wLJ^[{e + p)~^Ui "'j, showing the 
evolution of the gravitational wave energy density relative to 
an estimate of the square of the final fluid shear stresses. 
We take the fluid length scale Lf to be the integral scale ^f. 
Some oscillation about the constant curve is caused by long- 
wavelength sloshing of the fluid or the infrared behaviour of 
the gravitational wave power, discussed later, but the striking 
feature is the scalable linearity of the signal across a factor 
of three for i?*. Only fluid contributions to the gravitational 
wave power are included here. The early-times steep growth is 
best explained by the violent behaviour when the two shocks 
overlap. This phase is not well explained by our random ve¬ 
locity field model. 


phase transition is also of order the time it takes 

for bubbles of average separation to collide. 

Finally, there are also scales set by the background 
cosmology: the Hubble rate at the phase transition 
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(and the Hubble length), and the gravitational constant 
G. The Hubble rate, the gravitational constant and the 
critical temperature are related via the Friedmann 
equation. Our simulations are performed in a Minkowski 
background, as the duration of the transition is assumed 
to be comparable to the Hubble time. Therefore Tc and 
G can be chosen independently. The role of G is purely 
to set the scale of of the gravitational perturbations. As 
mentioned earlier, we use units Tc = 1 and G = 1. 

The observable of interest is the gravitational wave 
power spectrum, expressed as a fraction of the total den¬ 
sity (|32). It is clear from that formula that the relevant 
scales are the fluid flow length scale (set by the average 
bubble separation), the fluid flow lifetime, and the Hub¬ 
ble rate. The power spectrum is determined by the ratios 
of the fluid flow lifetime to the Hubble time, and the fluid 
flow scale to the Hubble scale. The role of the bubble wall 
width is to provide a short-distance cut-off on the power 
spectrum. A physical transition has i < Rc <Si R*, with 
of order 10“^ at the electroweak scale. Our sim¬ 
ulations assume that the bubble separation is much less 
than the Hubble length, and that the transition rate is 
much larger than the Hubble rate, so that expansion can 
be neglected. 

The fluid flow lifetime affects only the amplitude of the 
acoustically generated gravitational waves. In our simu¬ 
lations one sees that after the transition has completed, 
the power spectrum grows linearly with time while main¬ 
taining its shape. Hence, apart from a trivial scaling, the 
relevant parameter for the gravitational wave power spec¬ 
trum is the fluid flow scale, provided that the wall width 
and the critical bubble size are much less that the bubble 
separation. The effect of too large a ratio £/i?* is that 
there is insufficient dynamic range to observe the power 
law behaviour of the power spectrum; the effect of too 
large a ratio Rc/R* is that there is insufficient time for 
the fluid flow to approach its asymptotic self-similar pro¬ 
file, which results in too low a value for Ui. It also tends 
to obscure the power law behaviour. We have seen in our 
simulations that the ratio £/i?* needs to be of order 10“^ 
in order to reliably distinguish the power law. Given our 
computing resources, this means we are not able to deter¬ 
mine the shape of the power spectrum at wave numbers 
much less than the peak. 

In order to test the approach to physical ratios we 
should explore a scaling of the parameters which shrinks 
the ratios £/R^, and Rc/R* to zero. Such a scaling was 
given in Eq. (IMl) . Its only effect is to alter the width and 
surface tension of the bubble wall, and hence shrink the 
size of the critical bubble and the bubble wall width inde¬ 
pendently of the bubble separation. We carried out a sim¬ 
ulation scaled with r = 2 (so that the bubble wall was half 
the width) and parameters given in Table |H] correspond¬ 
ing to the deflagration, and compare the resulting grav¬ 
itational wave power spectra in Figure 1111 The power 
spectra are substantially similar, but the k~^ power law 
is clearer in the scaled run where i/R^, is smaller. This is 
consistent with our discussion above, and lends further 


ti=0.2 T , V =988, r=2500/r 

' c b c 



FIG. 11. Power spectra with Nh = 988 comparing the weak 
phase transition parameters and friction rj/Tc = 0.2 with the 
results from an equivalent run with the scaled parameters. 
For clarity, only the power spectra at the end of the phase 
transition {t — 2500/Tc) are shown. Note the y-axis scale 
is different to that used in Fig. [S] in order to highlight the 
differences between the two power spectra. 

confidence to our identification of the index of the power 
law in this case. 

Note that the scaling of Eq. (16411 also reduces the sur¬ 
face tension a as it reduces the bubble wall width, and 
hence the relative contribution of the scalar field to the 
total gravitational wave source tensor (l24l) . as the fol¬ 
lowing argument makes clear. The scalar field’s source 
tensor is proportional to the product of a with the 
area per unit volume of the phase boundary, and the 
area per unit volume is at most of order 1/i?*, which is 
unaffected by the scaling. Hence —>■ r~^T^. At the 

same time, the scale of the fluid source tensor is set by 
the latent heat of the transition £, which is independent 
of r. Hence the the relative importance of the scalar field 
to the fluid goes as as it decreases towards physical 
values. 

This is consistent with the argument given in the in¬ 
troduction that the scalar field contributes negligibly to 
the gravitational waves, as the ratio of the energy in the 
scalar field to the energy in the fluid goes as the ratio of 
the volume in the phase boundary to the total volume, 
which is i/R^:. 

This argument assumes that the bubble walls travel 
at constant speed, so that the effective surface energy is 
constant. Hence if the bubble walls are weakly coupled 
to the plasma, they can continue to accelerate until they 
collide [13 . In this “run-away” scenario, scalar fields can 
contribute importantly to the gravitational radiation. 


VI. DISCUSSION AND CONCLUSIONS 

In this paper we have reported on new numerical sim¬ 
ulations of the production of gravitational radiation at 
a first order phase transition in the early universe. Fol¬ 
lowing standard methods, we model the contents of the 
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universe as a scalar order parameter coupled to a rela¬ 
tivistic fluid, with a thermal effective potential © and 
dissipative coupling dH). This model captures the essen¬ 
tial physics of the transition, which proceeds by the nu- 
cleation and growth of bubbles of the low temperature 
phase. 

The most important parameters of the transition are 
the latent heat density relative to the total energy density 
ax, which characterises the strength of the transition, 
the bubble wall velocity Uw, which is determined by ax 
and the field-fluid coupling ry, and the bubble nucleation 
rate parameter (3, which determines the average bubble 
separation. 

Most of our simulations are carried out at ax — 0.01, 
the order of magnitude expected at an electroweak tran¬ 
sition, from which we can extrapolate to other values. 
We check our extrapolation with a smaller number of 
simulations at ax — 0.1, and with a scaling argument 
which changes parameters in the potential without af¬ 
fecting ax- We simulate for a range of phase boundary 
speeds Ww, covering deflagrations and detonations. In¬ 
stead of fixing the bubble nucleation rate parameter /3, 
we directly fix the average bubble separation ii* by nu¬ 
cleating iVb = V/R^ bubbles simultaneously. 

We concentrate on the gravitational waves generated 
by the fluid motion, as the vast majority of the latent 
heat of the transition is transformed into thermal and ki¬ 
netic energy of the fluid. We show that the gravitational 
wave density parameter (1431) is proportional to the fourth 
power of the mean square fluid velocity, the ratio of life¬ 
time of the source to the Hubble time, and the ratio of 
length scale of the source to the Hubble length. 

Our results confirm those of a more limited set of sim¬ 
ulations reported in Ref. [s^. The fluid kinetic energy is 
mostly in the form of sound waves generated by the com¬ 
pression or rarefaction of the fluid around the advancing 
phase boundary. Some rotational flow is generated by 
the collisions, but at a subdominant level. The sound 
waves remain for as long as we simulate, long after the 
phase transition completes. 

It was shown in Ref. [s^ that when viscosity is included 
the viscous damping time is much longer than the Hubble 
time for most phase transitions of interest. It was argued 
that the lifetime of the source, the shear stress generated 
by the sound waves, is approximately the Hubble time. 
In this paper we detail the calculation which shows that 
the lifetime parameter Tv, controlled by the decay and 
decorrelation of the shear stresses, is in fact exactly the 
Hubble time. 

The length scale of the source, app roximately the av¬ 
erage bubble separation in Ref. )34| . is here measured 
directly from the fluid flow. With this refinement, we 
show that the proportionality constant How in the grav¬ 
itational wave density parameter equation (1431) varies lit¬ 
tle between phase transitions with different strength and 
bubble wall speeds. Indeed, our measurements show that 
SttHgw = 0-8 ± 0.1, where the uncertainly is the root 
mean square fluctuation between simulations. 


Our new simulations are carried out on larger lattices, 
and give a wider dynamic range between the physical 
scales set by the average bubble separation and the bub¬ 
ble wall width. We further widen the dynamic range by 
nucleating all bubbles at the same time, at the slight 
cost of introducing “ringing” in the velocity power spec¬ 
trum. With the increased dynamic range we are able to 
establish clear power laws for both velocity and gravita¬ 
tional wave power spectra between the physical scales. 
For the transitions with = 0.44 or below, they are 
k~^ and k~^ respectively, where k is the wave number, 
and steeper for the transitions with higher bubble wall 
speeds. In order to discern these power laws, we show 
it is important that the fluid velocity profile around the 
advancing bubble wall has sufficient time to approach its 
asymptotic self-similar form. 

The k~^ (or steeper) power law for gravitational waves 
contrasts with the prediction of k~^ from the standard 
envelope approximation, which assumes that all the en¬ 
ergy in the system is concentrated in a thin shell at the 
bubble wall, and that the radiation is produced only 
when the shells interact. We see signs of a k~^ power 
spectrum generated by the scalar field in the initial phase 
of bubble collision, but this component is subdominant 
in our simulations, and would be completely negligible 
when extrapolated to the scale separation in a thermal 
phase transition. 

The envelope approximation generically predicts far 
less gravitational radiation than is actually produced. 
This under-prediction stems from the incorrect modelling 
of the source as being the colliding bubble walls. In¬ 
stead, the main source is the overlapping sound waves 
which are left behind after the transition has completed. 
We argued in [s^ that this means that the gravitational 
wave energy density is boosted by the ratio of the life¬ 
time parameter of the shear stress to the duration of the 
collision, which goes parametrically as . In 

this paper we studied the numerical factor in this ratio 
by a careful comparison of the quantities in the envelope 
approximation formula (1491) to the acoustic generation 
formula (l43l) . We show that the numerical factor is of 
order unity, and hence we can confirm that the gravita¬ 
tional wave signal is boosted by the ratio of the Hubble 
time to the phase transition duration, which is two orders 
of magnitude or more for a typical first order electroweak 
transition [40l. Id^. 

Our simulations shed new light on gravitational waves 
from phase transitions in the early universe. They show 
that the envelope approximation needs to be replaced, 
both as a model and as a formula. Instead, we should 
model the gravitational wave generation in terms of over¬ 
lapping sound waves. 

With this new “acoustic” model of gravitational wave 
generation, we have developed a quantitative under¬ 
standing of the gravitational wave density parameter 
(P)) . as a function of the mean square fluid velocity and 
the mean bubble separation. We can estimate the mean 
fluid velocity from hydrodynamic considerations , and 
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the mean bubble separation from the nucleation rate pa¬ 
rameter /3 [13 ■ We have numerically determined that 
the gravitational wave power spectrum is a power law 
on the high wavenumber side of the peak, and shown 
that it is steeper than the k~^ indicative of a vacuum 
transition. Hence potential future observations of such a 
gravitational wave spectrum will allow us to distinguish 
between a thermal and a vacuum transition. 

Much remains to be done. We noted that we need 
larger simulations to trace out the shape of the power 
spectrum at wavenumbers lower than the peak value, 
and to determine the index of the power spectrum for 
the transitions with faster bubble walls. They may also 
help in the search for bubble wall instabilities identified 
in [571 - 1^ . We also need to develop a theoretical under¬ 
standing of the shape or the power spectrum, and most 
importantly making accurate quantitative predictions for 
future gravitational wave observatories. 


Appendix A: Gravitational wave production in an 
expanding universe 

In this appendix we show how our discussion of the 
generation of gravitational waves in Section IIIII is modi¬ 
fied in an expanding radiation-dominated cosmology. We 
write the metric 


with = — 1. If we write 


C/O = W/a, f/‘ = WV^/a, (A3) 


then = 1/(1 — and resembles the Minkowski 
space 3-velocity. 

Indeed, it can be shown that the relativistic fluid 
equations in the radiation era are exactly the same as the 
Minkowski space equations, when the fluid variables are 
appropriately scaled. So writing 


4 4 

E = ^E, Z, = 

ai ai ’ 


(A4) 


with a* some reference scale factor, the equations for 
the tilde fields and are identical to the fluid parts 
of Eqs. (fTTll and (TT^ . with t interpreted as conformal 
time, and a;* as comoving coordinates. Hence our simula¬ 
tions need no adaptation for an expanding universe in the 
acoustic phase, where only fluid variables are active, and 
the expanding universe energy-momentum tensor can be 
obtained by multiplying by appropriate powers of the 
scale factor a. For example, the relevant quantity for 
gravitational wave generation is the fluid source tensor 
with both indices down (l24ll . This can be written 




(Al) 


with 77 representing conformal time. In the acoustic 
phase, after the scalar field has reached its ground state, 
the energy-momentum tensor takes the standard ideal 
fluid form 


= (e+p)t/^C/,+p5't, 


(A2) 


(k, 77 ) = 


02 ( 77 ) 


^il(k,?7), 


(A5) 


where we choose a* to be the scale factor at the phase 
transition time 77 *, and rb represents the source tensor 
obtained from the fluid evolution in scaled coordinates 
H*, E and Zi. 


We then see that the FLRW version of Eq. m is 


n^(fc, ^ 1 ,^ 2 ) — 


nttS- 


02(771)02(772) 


[{e+p)UfYLiIl^{kLf, kr]i,kr]2): 


(A 6 ) 
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where e, p, C/j and H are the values of the energy density pressure and mean square velocity in the scaled fields, Lf 

is the comoving length scale of the fluid perturbations, and II 2 is the same dimensionless function as in (IHTl) . We 
—2 

emphasise that e, p, C/f and H are those measured in our Minkowski space numerical simulations. 

In the metric (lAlll . the solution to the radiation era field equation for the tensor mode (|27l) is modified to 


/iy (k, 77) = (I67rG)Aij,fci(k) 


r'' sin[fc(77 - 77 ')] 0 ( 77 ') 

dr] -- ), 


k a{r]) 

and the definition of the gravitational wave energy density power spectrum (I23|) becomes 
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(A 8 ) 


dlogk 327rG 27r2 

Writing x = k{r]i -|- 772)72 and z = k{r]i — 772 ), and using the radiation era scale factor 0 ( 77 ) = (ji/rjYa*, the spectral 
density of h can be written 


Phik,v) = 
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where z± = ±2(a; —fcr;*) for x < k{r] + r]^,)/2 and z± = ±2{kr] — x) for x > k{r] + rjt.)/2. . Making the assumptions that 
the autocorrelation time of the fluid perturbations is small compared with the Hubble time, so that the integrand is 
negligible when z approaches x, and that the correlations are approximately stationary (independent of x) over the 
domain of integration we find 


Pl^(k,ri) ~ 


a\r]) V 




1 2 


167rG(e + p)C/f (a*7y*)(a*A: 


/ 

J — ( 




(AlO) 


We see that the expression has the same form as Eq. (1551) . but with the Minkowski time replaced by 0 * 77 * (the physical 
Hubble time) and the Minkowski wavenumber replaced by its physical value at time 77 *. We also see the correct redshift 
factor for the gravitational radiation. 

Hence we can immediately write down the analogue of Eq. (42) 


d^Gwjk) 

d\n{k) 


3{l + wfUi (iJ,Lf*) 


{kLtf 

277^ 


PcwikLf), 


(All) 


where Lf means the physical length scale at 77 ,, while PowikLf) remains its Minkowski space version (1371) . with k 
and Lf interpreted as comoving quantities. We therefore learn that the effective lifetime of the source for gravitational 
waves is precisely the Hubble time. 

Note that this effective lifetime is not the lifetime of the acoustic waves themselves, whose density perturbation 
continues to oscillate with constant amplitude in the absence of dissipation. Instead, it appears as a result of a 
combination of the expansion damping and decorrelation of the shear stress. To see the effect of the decorrelation, 
let us consider shear stress correlations behaving as nocos( 2 ;), in which case the z integrand in (IA10|) would always 
be positive, representing the largest possible growth rate for the gravitational wave power spectrum. In this extreme 
case, representing the effect of expansion damping alone, the factor (a* 77 *)(a*fc“^) in (lAlOl) would be replaced by 
as ( 0 * 77 *)^ ln^( 77 / 77 *). The decorrelation of the shear stresses cuts off the z integral, removing the logarithms and 
replacing one factor of the Hubble time with a factor of the wavelength. 


Appendix B: Simulations with a range of bnbble 
nucleation times 

In Section m we claim that the large oscillations in 
the power spectra presented in the main body of the pa¬ 
per are due to the bubbles being nucleated at the same 
time. There are many thin fluid shells of the same ra¬ 
dial size that contribute to the power spectrum, and the 
Fourier transforms of both the fluid velocity E* and un¬ 
projected metric perturbations hij will be in phase for all 
such shells. Here we demonstrate that such oscillations 
are damped when the bubbles are nucleated over a more 
widely spread period. 

We ran a single simulation with our “weak” poten¬ 
tial parameters and 77 /Tc = 0 . 2 , nucleating bubbles with 
a rate parameter P = 0.01. However, we capped the 
number of bubbles to approximately 1000 (in the end we 
nucleated 1002 ), leading to an abrupt cutoff in the expo¬ 
nential growth of the number of bubbles rather than the 
full double exponential seen in Refs. [13, The last 
bubble is nucleated shortly before t = 2QQ/Tc. 

While the number of bubbles and the spread of nucle- 
ation times will change the parameters we studied in the 
main body of the paper, they are sufficiently close for the 
purposes of this appendix. In Fig. [T2]we show the fluid 
velocity power spectrum (in the same manner as Fig. El) 
and demonstrate that the oscillations are considerably 
damped compared to the equivalent plot in Fig. El top 
right. 


ri=0.2 T^, A^|^=1002, unequal nucleation time 



FIG. 12. Velocity power spectrum with a range of nucleation 
times. Compared to the similar simulation results in Fig. El 
top right, the oscillations in the fluid power are significantly 
more damped. This demonstrates that the strong oscillatory 
behaviour in the fluid velocity power spectra is due to the 
equal nucleation time. 


Appendix C: Transverse modes of the velocity field 
are negligible 

Here, we show that the power in the transverse modes 
of the fluid velocity is significantly smaller than that in 
the longitudinal modes, supporting our claim that the 
fluid perturbations are best characterised as an essen¬ 
tially linear superposition of sound waves. 
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FIG. 13. Comparison of longitudinal (solid) and transverse 
(dashed) fluid velocity power spectra for = 988, r\lTc. = 
0.2. The time intervals are 500/Tc as in Fig. [3 (c.f. top 
right figure). The transverse part of the fluid velocity power 
is many orders of magnitude smaller than the longitudinal 
part and does not play a significant role in gravitational wave 
production. Note the extended scale compared to Fig. [7] 

One way of quantifying this is to study the RMS fluid 
velocity when only transverse motion is taken into ac¬ 
count. We quantify this with C/f TablellTl It can be 

seen that the transverse modes contribute at most 5—10% 
to the RMS fluid velocity I/f. This ratio is greatest for 
the simulations with the largest iVb (and hence smallest 
R*), so we would expect in a realistic scenario (with i?* 
considerably larger) that the transverse modes would be 
very small indeed. 

The maximum value oiUf (which is the value quoted 
in Table HI]) occurs at approximately the same time as the 
maximum value of I/f, at around the conclusion of the 
phase transition before slowly decreasing due to discreti¬ 
sation effects. This behaviour can be seen in Figure [SJ 

Finally, we show in Figure[T3]the transverse fluid veloc¬ 
ity power spectrum (dashed lines) for a typical simulation 
alongside the longitudinal power (solid lines). It can be 
seen that it is significantly smaller than the longitudinal 


power (a ratio of around 10 ^). 

It has been predicted that turbulent fluid motion would 
develop after a phase transition such as the one under 
study in this paper. Fluid turbulence is generally stud¬ 
ied in incompressible (i.e. rotational) flows, with energy 
transport from a larger forcing scale to a smaller dissi¬ 
pation scale, leading to the formation of characteristic 
power laws. No clear power law can be seen in the trans¬ 
verse velocity power spectrum, which leads us to believe 
that turbulence is not a feature of our simulations. 

It has recently been suggested that turbulence can de¬ 
velop in the acoustic perturbations , giving rise to an 
inverse cascade (transfer of power to longer scales). We 
see no signs of the length scale of the acoustic oscillations 
changing once the transition is complete, and the velocity 
power spectrum does not change its form significantly. 

While it is possible that turbulence develops at larger 
Reynolds numbers than we have access tc0, it is clear 
that it has no significance for gravitational radiation in 
the relatively weak transitions we study. 
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